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Abstract. In the resonance model, high-frequency quasi-periodic oscillations (QPOs) are supposed to be a consequence of 
nonlinear resonance between modes of oscillations occurring within the innermost parts of an accretion disk. Several models 
with a prescribed mode-mode interaction were proposed in order to explain the characteristic properties of the resonance in 
QPO sources. In this paper, we examine nonlinear oscillations of a system having two degrees of freedom and we show that 
this case could be particularly relevant for QPOs. We present a very convenient way how to study autoparametric resonances 
of a fully general system using the method of multiple scales. We concentrate to conservative systems and discuss their 
behavior near the 3:2 parametric resonance. 
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1. Introduction 

In the resonance model l|AbnmTO\¥icz& Kluzniakl l200lt 
iKluzniak & Abramowic jl2000t lKatoll2003l) ~ there is a natu- 
ral and attractive possibility of explaining the observed ratio- 
nal ratios of high-frequency QPOs as a consequence of non- 
linear coupling between different modes of accretion disk os- 
cillations. The idea has been pursued in several p apers (re- 
cently, e.g. lAbramowicz et all2003URebuscol2004 . 

Specific models invoke particular physical mechanisms. 
Some models can be almost immediately comprehended 
as distinct realizations of the general approach discussed 
here - for example, various formulations of the orbiting 
spot model (Schnittman & Bertchinger 2004) or the models, 
where QPOs are produced by the magn etically d riven res- 
onance in a diamagnetic accretion disk (ILaill 19991) - while 
other seem to be more distant from the view presen ted herein 
- e.g. the transition layer model dTitarchukl l2002). an inter- 
esting idea of p-mod e oscillations of a small accretion torus 
dRezzola et alJl2003h or the model of blobs i n an accretion 
disc (see e.g. lKarasll999tLi & Naravan|2004 a nd references 
cited therein). Also in this context, Kato ( 2004) discussed the 
resonant interaction between waves propagating in a warped 
disk, including their rigorous mathematical description. In- 
stead of pursuing a specific model, here we keep the discus- 
sion as general as possible, aiming to implement the formal- 
ism of multiple scales. Indeed, we show that there is unques- 
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tionable appeal in this approach which offers some additional 
insight into generic properties of resonant oscillations. 

Some properties of an accretion disk oscillations can be 
discussed within the epicyclic approximation of a test particle 
on a circular orbit near equatorial plane. Suppose that angular 
momentum of the particle is fixed to a value £. The effective 
potential Ug(r, 9) has a minimum at radius ro, corresponding 
to the location of the stable circular orbit. An observer mov- 
ing along this orbit measures radial, vertical and azimuthal 
epicyclic oscillations of a particle nearby. Since the angular 
momentum of the particle is conserved, only two of them - 
radial and vertical - are independent. The epicyclic frequen- 
cies can be derived from the geodesic equations expanded to 
the linear order in deviations Sr = r — ro and 89 = 8 — ir/2 
from the circular orbit. We get two independent second-order 
differential equations describing two uncoupled oscillators 
with frequencies oj r and ojg, which are given by the second 
derivatives of effective potential Ut(r, 9). 

In Newtonian theory, uj r and uig are equal to the Keple- 
rian orbital frequency VLk- This is in tune with the fact that 
orbits of particles are planar and closed curves. The degen- 
eracy between two epicyclic frequencies can be seen as a re- 
sult of scale-freedom of th e Newtonian gravitational potential 
jAbramowicz & Kluzniakl2003l) . In Schwarzschild geometry 
this freedom is broken by introducing the gravitational ra- 
dius r g = 2GM/c 2 . The degeneracy between the vertical 
epicyclic and the orbital frequencies is related to spherical 
symmetry of the gravitational potential, which assures the 
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existence of planar trajectories of particles. All three frequen- 
cies are different in the vicinity of a rotating Kerr black hole. 

In addition, when nonlinear terms of geodesic equations 
are included, the two oscillations in r and 9 directions be- 
come coupled and variety of new phenomena connected 
to nonlinear nature of the equations appear. This rich phe- 
nomenology includes frequency shift of observed frequen- 
cies with respect to eigenfrequencies, presence of higher har- 
monics and subharmonics, drifts and parametric resonance. 
The first three are connected to nonlinear oscillations of each 
mode and the last one comes from the coupling between two 
modes. 

2. Expansion via multiple scales 

We study nonlinear oscillations of the system having two de- 
grees of freedom, i.e., the coordinate perturbations 5r and 
59. The oscillations are described by two coupled differen- 
tial equations of the very general form 

Sr + u% 5r = J 1 f r {8r, 59, 5r, 59), 



+ u 2 56 = ujj f g (5r,59,5r, 



(1) 
(2) 



Suppose that the functions f r and fg are nonlinear, i.e., their 
Taylor expansions start in the second order. Our another as- 
sumption is that these functions are invariant under reflection 
of time (i.e., the Taylor expansion does not contain odd pow- 
ers of time derivatives of Sr and 59). As we shell see later, 
this is related to the conservation of the total energy in the 
system. Many authors stud ied such systems with a particular 
form of functions / and g (jNavfeh & M ookl 19791) . however, 
in this paper we keep the discussion fully general. 

We seek the solutions of the govern ing equations in 
the form of the multiple-scales expansions (Navfeh&Mook 
Il979h 



5r(t,e) 



4 

£ 



e n r n (T^), 50{t,e) = Y J t n MT t * 



(3) 



where several time scales T M are introduced instead of the 



physical time t, 

T^ = eH, M = 0,l,2,3. 



(4) 



The time scales are treated as independent. It follows that in- 
stead of the single time derivative we have an expansion of 
partial derivatives with respect to the Ij, 

j t = D + eD 1 + e 2 D 2 + e 3 D 3 + 0(e 4 ), (5) 
d 2 

— = D 2 + 2eD D 1 + e 2 (D 2 + 2D Q D 2 ) + 

2e 3 {D D 3 + D 1 D 2 ) + O(e 4 ), (6) 

where = d/dT^. 

We expand the nonlinear functions f r and fg into the Tay- 
lor series and then we substitute the expansions (|3jl, l|5} and 
(|6}. Finally, we compare the coefficients of the same powers 
of e on both sides in the resulting couple of equations. This 
way we get a set of linear second-order differential equations 
that can be solved successively - the lower-order terms of the 



expansion (0 appear as forcing terms on the right-hand sides 
of the equations for the higher order approximations. 

In the first order we obtain equations corresponding to the 
linear approximation 

(D 2 +oj 2 ) ri =0, (Dl+wDe^O. (7) 
with the solutions 

x x = A r (T x ,T 2 , T 3 )e l ^To + CCj (8) 
9 1 =A g {T 1 ,T 2 ,T 3 )e l ^ T °+cc. (9) 

The complex amplitudes A r and Ag generally depend on the 
higher time-scales. 

The solutions l|8} and l|9} substituted into the quadratic 
terms on the right-hand sides of the second-order differential 
equations produce terms that oscillates with frequencies 2ui r , 
2u)g and uig ± ui r . When the frequency ratio u r /wg is far from 
1:2 and 2:1 the particular solutions r 2 and 9 2 describe higher 
harmonics to the linear-order oscillations r\ and 9\ . Hence, 
the presence of higher harmonics in the power-spectra is a 
signature of nonlinear oscillations. Their frequencies and rel- 
ative strengths with respect to the main oscillations provide 
us an usefull informations about nonlinearities in the system. 

In addition, the right hand sides of the second order equa- 
tions contain terms proportional to e luJrT ° and e tuJeTo that os- 
cillates with the same frequency as the eigenfrequency of 
the oscillators. These terms produces secular grow of the 
amplitudes of the second-order approximations r 2 and 9 2 
and causes nonuniform expansions (|3jl. Eliminating them we 
get the solvability conditions for the complex amplitudes 
A r {T x ,T 2 , T 3 ) and A e (T 1 ,T 2 ,T 3 ) t hat give us the evolutio n 
of the system on longer time-scales JNavfeh & Mookl 19791) . 

When the eigenfrequencies are in 1:2 or 2:1 ratio we ob- 
serve qualitatively different behavior related to the autopara- 
metric resonance. In that case the right hand sides contains 
additional secular terms and the solvability conditions take 
different form. Different resonances occur in different orders 
of approximation. The possible resonances in the third order 
are 1:3, 1:1 and 3:1 and 1:4, 3:2, 2:3 and 4:1 in the fourth or- 
der 1 However, if the governing equations remain unchanged 
under the transformation 59 — » —59 (i.e., the system is re- 
flection symmetric) the only autoparametric resonanc es tha t 
exists in the system are 1:2, 1:1, 1:4 and 3:2 (Rebusco 2004) 

3. The 3:2 autoparametric resonance 

Let us consider oscillations of a conservative system eigen- 
frequencies of which are close to 3:2. The time behavior of 
the observed frequencies ui* and ojg and amplitudes a r and 
ag of the oscillations can be found from the solvability con- 
ditions imposed on the complex amp litudes A r (Ti,T 2 ,T 3 ) 
and Ag(Ti . T,.n) JHorak et al.l2004 



D x A r = D 2 A g = 0, 

D 2 A r = [ Kr \A r \ 2 + Kg\Ag\ 2 ] A r , 

D 2 Ag = ^[\ r \A r \ 2 + \g\Ag\ 2 ] Ag, 
1 The ratio n : m refers to the eigenfrequency ratio cue : ui r 



(10) 

(11) 

(12) 



3.2 Integrals of motion 
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D 3 A r = - l -u r a{AlyA*e- i ^ T *+°* T *\ 



3^9 



-u! 9 f3A 3 r A* g e 



i(a 2 T 2 +<T 3 T 3 ] 



(13) 



(14) 



In the fourth order we eliminate also terms which become 
secular in when 3oj r « 2uog. We describe vicinity of the res- 
onance by the detuning parameters 02 and 03 introduced ac- 
cording to 

3uj r = 2uj e + e 2 <7 2 + £ 3 ct 3 . (15) 
The term e<j\ is missing, because the complex amplitudes de- 
pends only on the second and the third time-scales. The solv- 
ability conditions describe evolution of the system in the most 
general way: the real parameters a, (3, K r , Kg, X r and Xg are 
given by the coefficients of the Taylor-expanded nonlinear 
functions f r and fg. 

Since A r and Ag are complex, the conditions dlOi — J14I 
represents 12 real equations. However few of them are trivial. 



By substituting the polar forms eA r 

pi<t>0 



■ e^ and eAg 



^age" t ' e , separating real and imaginary parts and introducing 
the unique time t the number of the equations can be reduced 
to four, 



a„ 



a 



auj. 

~l6 Up 

(3ujg 



a 2 a 2 sin 7, 



16 
~2 

UJg 



a" a e sin 7, 



= ~— \n r a 2 + Kg a 2 ] 



Xg a g 



auj r 
TJT 

P^g 

16 



a p a B cos 7, 



cos 7, 



(16) 
(17) 
(18) 
(19) 



where we introduced the phase function ^(t) = 2(f)g(t) — 
3(f> r (t) — at and the unique detuning parameter a = e 2 (T2 + 
e 3 er 3 . The equations fl!6l > and d!7t describes the slow evolu- 
tion of the amplitudes of oscillations and additional long-term 
behavior of the oscillation phases is given by equations (II 81 
and ( 1191 . These equations give us the frequency-shift of the 
observed frequencies uo* and oj* with respect to the eigenfre- 
quencies uj r and ujg, respectively, 

UJ* = UJ r + <f) r , UJg = UJg + (j)g. (20) 

The two equations Jl 8I > and dl9l can be replaced by a sin- 
gle differential equation for the phase function, 



7 



UJg 
"4 



\i, r a r 



\iga 2 g 



y 



f3a 2 ) cos 7 



,(21) 



were we used the fact that near the resonance u r w (2/3)u>g 
and we defined \i r — k t — X r and fig — Kg — Xg. 

3.1. Steady-state solutions 

Steady-state solutions are characterized by constant ampli- 
tudes and frequencies of oscillations. Such solutions repre- 
sent singular points of the system governed by equations dl6> , 
(O and (ED. 

It is obvious from equations fl!6l > and dl7t that the con- 
dition a r = ag = can be satisfied (with nonzero ampli- 
tudes) only if sin 7 = (identically at all times), and thus 
also 7 = 0. In that case equation i2l\ transforms to the alge- 
braic equation 
a 1 



fi r a r + figa'f 



Or 

2 



aag 



(3a 2 



(22) 
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Fig. 1. Comparison between an analytical constraint 1271 
and the corresponding numerical solution of the system 
lAbramowicz et alJ 112003 ). Each point corresponds to the am- 
plitudes of the oscillations at a particular time. On the other 
hand, from the discussion of equation J27l > we know that these 
points must lay on an ellipse, whose shape is determined by 
the multiple-scales method. 



The left-hand side can be expressed using the eigenfrequency 
ratio R = ojg /ui r as 



UJg R\ 2 

Then we get 



(23) 



R 



3 3 3 

2 ~ 16 ^ r °^ + M<?a ^ ± 32 ar ~ ^""^ ' ^ 



were we neglected terms of the fourth order. Note that the 
lowest correction to eigenfrequencies is of the second order - 
for given amplitudes a r , ag the steady-state oscillations occur 
when the eigenfrequency ratio departs from 3/2 by deviation 
of order of a 2 . 

The relation between observed frequencies of oscillations 
oj*, uj* and eigenfrequencies oj r , ujg are given by the time 
derivative of phases <j) r and cf>g. We can find simple relation 
between observed frequencies and the phase function 



Suit — 2UJ* 



3uj r — 2ujg 
a + (3(f) r - 



+ (30,. 
2fo) -- 



2(fjg) 
-7- 



(25) 



For steady state solutions 7 = 0, and thus observed frequen- 
cies are adjusted to exact 3:2 ratio even if eigenfrequencies 
depart from it. 

3.2. Integrals of motion 

The method of investigation of time-dependent behavior of 
the system is analo gical to the cas e of 1:2 resonance as ex- 
amined bv lNavfeh & Moot] \ 19791) . The oscillations are de- 
scribed by three variables a r (t), ag(t) and j(t) and three 
first-order differential equations (I16> . dl7> and (121 1 . How- 
ever, the number of differential equations can be reduced to 
one because it is possible two find two integrals of motion. 
Our discussion will be 
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3.3 Analytical results 



Consider equations (1161 and d!7i . Eliminating sin 7 from 
both equations we find 

d 

eft" 

and thus 



-{al + voe) = 



(26) 



a^. + vag — const = E, 



(27) 



(28) 



where we defined 

auj r 2a 

When v > 0, the both amplitudes of oscillations are bounded. 
The curve [a r (t), ag(t)] is a segment of an ellipse. The con- 
stant E is proportional to the energy of the system. On the 
other hand, when v < 0, one amplitude of oscillations can 
grow without bounds while the second amplitude vanishes. 
This case correspon ds to the presence of an regenerative ele- 
ment in the system (Navfeh & Mook 1979). The correspond- 
ing curve in the (a r , ag) plane is a hyperbola. In further dis- 
cussion we assume that v > 0. 

In order to verify that the the energy of the system 
is conserved, we numerically integrated governing equa- 
tion and J3 fo r the o ne particular system discussed by 
Abramowi cz et alJ d2003h . The comparison is in Figure 
The numerical and analytical results are in a very good agree- 
ment. 

The second integral of motion is found in following way. 
Let us multiply equation i2l\ by ag . Then we obtain 



agj 



Ulg 2 LOg , Log ^> 

-crag + —fj, r a r ag + —^ga g + —aa r a g cos 7 
4 4 8 



LUt > a 3 

—pa r ag cos 7. 



(29) 



Changing the independent variable from t to ag and multiply- 
ing the whole equation by dag we find 



8(7 



a r a g d(cos 7) + — — d(a g ) —a r agd(ag) 7T^( a 



ftujg 



ft 



ft 



2a 



— a~ ara 6 cos 7^ a # + 2a^ag cos 7dag = 0. (30) 

P 

The equation J27I implies 

a r da r 

agdag = . (31) 

v 

With the aid of this relation equation ( I30i takes the form 

3a^ag cos jda r + 2a^ag cos jdag + a^agd(cos + 

+^d(a 2 e ) + ^d(ai) - f^d(4) = 0. (32) 

ftiUg ftv (j 

The first three terms express the total differential of function 
— a^ag cos 7. Hence, the above equation can be arranged to 
the form 

d ( ala 2 9 cos 7 + i^ag + j^ a A _ A = (33) 



ftuog a ftv r (3 



In other words, 



» 2,^4 V0 4 



a r a g cos 7 + — — a g + —a — a g = const = L (34) 

ftijjg ftv 

is another integral of equations JT6l . flTl and (I2TV 




Fig. 2. Example of the (7,£)-plane for the system close to 
the 3:2 resonance. The oscillations are coupled by nonlinear 
functions f p and fg [see equations Q and Q]. These func- 
tions give us values of the constants a, ft, K r , Kg. \ r and \g. 
The thick solid line is the separatrix dividing the librating and 
circulating trajectories The blue dotted line connects points 
where 7 = 0. The example is for values a = [3 = n r = \g = 
l, K g = Xg = 2,S = 0.1 and a = -0.165. 



3.3. Analytical results 

Knowing two integrals of motion, we are able to find one 
differential equation which governs the time-evolution of the 
system. 

First, the amplitudes a r and ag are not independent be- 
cause they are related by equation J27i . To satisfy this rela- 
tion, let us define new variable £(t) by 



= £E, al = {!-£)-. 

v 



(35) 



For present moment, we ignore the time dependence by 
considering projections of solutions into the (7. £) -plane. For 
a fixed energy E of oscillations, the system follows curves 
of constant L. Hence, the trajectories in the (7, £)-plane are 
given by equation 

L(j. £) = const. (36) 

An example of the phase-plane is given in Figure [2] There 
are two types of trajectories [£(t),j(t)]: the circulating tra- 
jectories take the full range < -f(t) < 2ir and the li- 
brating trajectories that are confined in the smaller range 
7i < j(t) < 72- The turning points on the librating trajec- 
tories correspond to 7 = 71 and 7 = 72- This division has 
an interesting consequences with respect to the observed fre- 
quencies of resonant oscillations. According to the relation 
( I25> the observed frequencies are in exact 3:2 ratio when the 
system pass through these points. On the other hand, the cir- 
culating trajectories do not contain any turning points and the 
ratio of observed frequencies is always above or bellow 3:2. 

The equation describing the evolution of £(t) can be de- 
rived in the following way. Let us multiply equation d 1 6t > by 
2a r and integrate it. We obtain 

d(4) 



dt 



" 3 2' 

—aj r a r ag sin 7 



(37) 



3.3 Analytical results 
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Fig. 3. The functions ±F{£) = ±(1 - £)£ 3/2 and the 
quadratic function G(£) whose second powers are first and 
second terms on the right-hand side of equation J40I . The be- 
havior of the system corresponds to £ in the interval £2] 
(denoted by the two dotted vertical lines) where the condition 
\F(Cl\ > \G(Q\ is satisfied. 



Then we express a 2 using £, and square it. We find 

8E N 2 



( 3 2 • \ ^ 

I a r a e sin 7 j . 



(38) 



The right-hand side of this equation can be expressed using 
equation J34i as 

/ 3 2 • \2 / 5 o\2 

(a r a e SI117J = [a r a e j — 

2 



8a 



L - - — a H - — a* + — a„ ) . (39) 



After the substitution into equation J38i and using the re- 
lations J35I . we get 



-(—) ? 



2 e 3 



(!-&'( 

fJ, r E 2 



8aE 

flvujg 

, A^ 2 



(i-O 

2 (1- 



2 ] 2 . 



(40) 



The equation of motion has a form 

/c 2^2 =F 2 (c) _ G 2 (e)i (41) 

where the /C 2 is a positive constant, = (1 — £)£ 3//2 and 
G(£) is a quadratic function coefficients of which depend on 
initial conditions through E and L. The motion occurs only 
for £ that satisfy > |G(£)|. The turning points, where 

£ changes its signature, are determined by the condition 

\F(0\ = \G(0\- (42) 

The functions ±F(£) and G(£) are plotted in Figure [3] 
Generally, the function G intersects the functions ±F in two 
points that corresponds to oscillating between the two 
bounds £1 and £2 given by condition (I42i . The radial and 
vertical mode of oscillations periodically exchanges the en- 
ergy. The amount of exchanged energy is given by AE/E — 
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Fig. 4. Time evolution of the amplitudes (top panel) and the 
observed frequencies, v* e = o;g/(27r) (middle panel) and 
v* = oj*/(2tt) (bottom panel). All quantities are rescaled 
with respect to the higher eigenfrequency vg . Amplitudes of 
the oscillations are anticorrelated because the energy is con- 
served. Observed frequencies are correlated because the sys- 
tem is in parametric resonance. 



£2 — £i- For some particular values of L and E only one in- 
tersection of ±F and G exists (the function G(£) touch one 
of the functions ±F(£))- the oscillations of the system cor- 
respond to the steady-state solutions discussed above. 

The period of the energy exchange can be find by integra- 
tion of equation iAQl 



T = iii?- 3 / 2 - 



This integral can be roughly approximated as 
16tt 



T 



(3ug' 



_ E -3/2_ 



(43) 



(44) 



However, near the steady state the period becomes much 
longer. 

The observed frequencies lo* and ojg are given by rela- 
tions J25b . They depend on squared amplitudes a 2 and a 2 ,. 
Since both a 2 and a 2 , depend linearly on also observed 
frequencies are linear functions of £ and are linearly corre- 
lated. The slope of this correlation ojg = Kuo* + Q is in- 
dependent of the energy of oscillations and is given only by 
parameters of the system, 



K 



cue \ r v - X g 

UJ r K r V — Kg 



(45) 



The slope of the correlation differs from 3:2, however the ob- 
served frequencies are still close to it. 
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3.4. Numerical results 

The equations d!6t . d!7> and fl21i were solved numerically 
using the fifth-order Runge-Kutta method with an adaptive 
step size. One of the solutions is shown in Figure |4] The 
top panel of the figure shows the time behavior of the ampli- 
tudes of the resonant oscillations. Since energy of the system 
is constant, amplitudes are anticorrelated and the two modes 
are continuously exchanging energy between each other. The 
middle and the bottom panels show the two observed frequen- 
cies that are mutually correlated. They are also correlated to 
one of the amplitudes. The frequency ratio varies with time 
and it differs from the exact 3:2 ratio, however, it always re- 
mains very close to it. Our numerical solution is in agreement 
with the general results obtained analytically in the previous 
section. 

4. Conclusions 

Although this paper was originally motivated by observations 
and models connected to high-frequency QPOs, our results 
are very general and can be applied to any system with gov- 
erning equations of the form Q and 0. 

The main result of the calculations is our prediction of 
low-frequency modulation of the amplitudes and frequencies 
of oscillations. The characteristic timescale of the modulation 
is approximately given by equation J44i . 

Because of the generality of our approach this fact have 
an interesting consequences in the context of QPO nature 
of which are unknown. Our result can be summarized in 
the following way: If the two quasiperiodic oscillations ob- 
served close to 3:2 ratio are produced by the autoparamet- 
ric resonance the frequencies and amplitudes of oscillations 
should be periodically modulated. This modulation appears 
as a separate peak at the modulation frequency and as side- 
ba nds to the main (li near) oscillation. In a separate paper 
by Ho rak et alJ J2004I) we pointed to possible connection of 
this modulation with the 'normal branch oscillations' (NBOs) 
that are often present together with QPOs. Specifically, we 
suggest that the correlation between the higher frequency and 
the lower amplitude, evident in Figure |4] is the same as was 
recently reported in Sco X-l bv lYu et alJll200ll) . We note that 
similar behavior was recently observe d also in the g alactic 
black-hole candidate XTE J1550-564 dYu et all2002l) . 
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